
setClass("AEP.frame", representation(
  data = "data.frame",
  pobs_data = "matrix",
  nb_obs = "numeric",
  nb_var = "numeric",
  EmpiricalCDF = "list",
  ECDF = "function",
  N = "numeric",  
  binary_matrix = "array",
  sum_ones = "array",
  Hmultiplier = "array",
  alpha = "numeric",
  iteration = "numeric",
  m = "array",
  bandwidth = "numeric",
  data.weighted_bandwidth = "matrix"
)
)

setMethod("initialize", "AEP.frame", function(.Object, data, alpha, iteration=1, bandwidth) {
  
  if (missing(data)) stop("syntax error")
  if (is.data.frame(data) == FALSE) stop("the parameter 'data' must be of type 'data.frame'") 	  
  
  if (missing(alpha)) alpha <- 2 / (dim(data)[2] + 1)
  else if (is.numeric(alpha) == FALSE) stop("the parameter 'alpha' must be of type 'numeric'")      
  if (alpha < 1 / dim(data)[2] || alpha >= 1) stop("the parameter 'alpha' must be in [1/d, 1)")   		
  
  if (is.numeric(iteration) == FALSE) stop("the parameter 'iteration' must be of type 'numeric'") 	 
  if (iteration < 1) stop("the parameter 'iteration' must be >= 1")			
  
  if (missing(bandwidth)) bandwidth = rep(0, c(dim(data)[2]))
  if (is.numeric(bandwidth) == FALSE) stop("the parameter 'bandwidth' must be of type 'numeric'")    
  if (length(bandwidth) != dim(data)[2]) stop("the parameter 'bandwidth' must be of same length as the number of variables in 'data'")    
  
  .Object@data <- data
  .Object@nb_obs <- dim(data)[1]
  .Object@nb_var <- dim(data)[2]  
  .Object@alpha <- alpha	
  .Object@iteration <- iteration
  .Object@N <- 2^.Object@nb_var
  .Object@binary_matrix <- array(0, c(.Object@N, .Object@nb_var))			
  .Object@sum_ones <- array(0, c(.Object@N))
  .Object@Hmultiplier <- array(0, c(.Object@N))
  .Object@bandwidth <- bandwidth  
  .Object@data.weighted_bandwidth <- t(.Object@bandwidth * t(.Object@data))  
  
  #Compute the ECDF of each variable
  EmpiricalCDF <- NULL	
  for (i in 1:.Object@nb_var) .Object@EmpiricalCDF <- c(.Object@EmpiricalCDF, ecdf(.Object@data[,i]))
  .Object@ECDF <- function(i, x) .Object@EmpiricalCDF[[i]](x)
  
  # Compute the pobs of data
  .Object@pobs_data <- pobs(.Object@data)
  #.Object@pobs_data <- array(0, c(.Object@nb_obs, .Object@nb_var))
  #for (i in 1:.Object@nb_obs) {
  #  for (j in 1:.Object@nb_var) {
  #    .Object@pobs_data[i,j] <- .Object@ECDF(j, .Object@data[i,j]) 
  #  }      
  #}
  #.Object@pobs_data <- as.matrix(.Object@pobs_data)
  
  #Initialize the matrix of binary and the sum of ones in each line
  for (i in 2:.Object@N) {
    .Object@binary_matrix[i,] <- .Object@binary_matrix[i-1,]
    j <- 1
    .Object@binary_matrix[i, j] <- .Object@binary_matrix[i, j] + 1
    while (.Object@binary_matrix[i, j] > 1) {
      .Object@binary_matrix[i, j] <- 0
      j <- j + 1
      .Object@binary_matrix[i, j] <- .Object@binary_matrix[i, j] + 1
    }
    .Object@sum_ones[i] <- sum(.Object@binary_matrix[i,])
  }
  
  #Initialize the multiplier of H()
  for (i in 1:.Object@N) {
    if ((.Object@nb_var  -  .Object@sum_ones[i]) %% 2 == 0) .Object@Hmultiplier[i] <- 1
    else  .Object@Hmultiplier[i] <- -1
  }
  
  # Initialize the vector m
  .Object@m <- array(0, c(.Object@N))
  for (j in 2:.Object@N) {
    if (.Object@sum_ones[j] < 1 / .Object@alpha) {
      if((1 + .Object@sum_ones[j]) %% 2 == 0) .Object@m[j] <- 1
      else .Object@m[j] <- -1
    }
    else if (.Object@sum_ones[j] == 1 / .Object@alpha) .Object@m[j] <- 0
    else {
      if((.Object@nb_var + 1 - .Object@sum_ones[j]) %% 2 == 0) .Object@m[j] <- 1
      else .Object@m[j] <- -1
    }
  }
  .Object
})

setMethod(f = "show", signature = "AEP.frame", definition = function(object) {
  cat("\n@nb_obs = number of observations = ", object@nb_obs, "\n\n")
  cat("@nb_var = number of variables = ", object@nb_var, "\n\n")
  cat("@data = [not displayed]\n\n")
  #print(object@data)
  cat("@pobs_data (pobs(data)) = [not displayed]\n\n")
  #print(object@pobs_data)
  cat("@ECDF(i, x) = empirical cumulative distribution function of var 'i' at 'x', i.e. P[var_i <= x]\n\n")
  cat("@N = pow(2, nb_var) = ", object@N, "\n\n")
  cat("@binary_matrix = \n")
  print(object@binary_matrix)
  cat("\n@sum_ones = sum of each row of the binary matrix = ", object@sum_ones, "\n\n")
  cat("@Hmultiplier = (-1)^(nb_var - sum_ones) = ", object@Hmultiplier, "\n\n")
  cat("@alpha = ", object@alpha, "\n\n")
  cat("@iteration = ", object@iteration, "\n\n")
  cat("@m =", object@m, "\n\n") 
  cat("@bandwidth =", object@bandwidth, "\n")
})

as.aep.frame <- function(data, alpha, iteration=1, bandwidth) {
  if (missing(data)) stop("syntax error")
  if (is.data.frame(data) == FALSE) stop("the parameter 'data' must be of type 'data.frame'")       
  if (missing(alpha)) alpha <- 2 / (dim(data)[2] + 1)
  else if (is.numeric(alpha) == FALSE) stop("the parameter 'alpha' must be of type 'numeric'")      
  if (alpha < 1 / dim(data)[2] || alpha >= 1) stop("the parameter 'alpha' must be in [1/d, 1)")     	
  if (is.numeric(iteration) == FALSE) stop("the parameter 'iteration' must be of type 'numeric'") 	 
  if (iteration < 1) stop("the parameter 'iteration' must be >= 1")	    
  if (missing(bandwidth)) bandwidth = rep(0, c(dim(data)[2]))
  if (is.numeric(bandwidth) == FALSE) stop("the parameter 'bandwidth' must be of type 'numeric'")    
  if (length(bandwidth) != dim(data)[2]) stop("the parameter 'bandwidth' must be of same length as the number of variables in 'data'")      
  new("AEP.frame", data, alpha, iteration, bandwidth)
}

VHmeasure <-function(frame, hc) {
  print("entering VHmeasure: ")
  if (missing(frame) || missing(hc)) stop("syntax error")
  if (class(frame) != "AEP.frame") stop("the parameter 'frame' must be of type 'AEP.frame'")
  if (class(hc) != "AEP.HCSimplex") stop("the parameter 'hc' must be of type 'AEP.HCSimplex'")
  if (length(hc@b) != frame@nb_var) stop("hc@b is too long")    
  Kn <- array(0, c(frame@nb_obs))  
  # variable to increase the execution speed
  hc.sum <- hc@b + hc@h  
  u.upper <- rep(1, frame@nb_var)
  u.lower <- rep(-1, frame@nb_var)
  total.prod <- prod(u.upper - 1/3 * u.upper^3 - u.lower + 1/3 * u.lower^3)
  for (i in 1:frame@nb_obs) {    
    u.upper <- (hc.sum - frame@data[i,]) / frame@data.weighted_bandwidth[i,]
    u.lower <- (hc@b - frame@data[i,]) / frame@data.weighted_bandwidth[i,]    
# Integrate version of the multivariate Epanechnikov kernel density function (multiplicative)              
    if (any(u.lower > 1, u.upper < -1)) Kn[i] <- 0    
    else {
      u.upper[u.upper > 1] <- 1
      u.lower[u.lower < -1] <- -1            
      if (all(u.upper == 1 & u.lower == -1)) Kn[i] <- total.prod
      else Kn[i] <- prod(u.upper - 1/3 * u.upper^3 - u.lower + 1/3 * u.lower^3)    
    }  
  }
  print("exiting VHmeasure")
  (1/frame@nb_obs) * (3/4)^(frame@nb_var) * sum(Kn) 
}


resolveAEP <- function(frame, vsum) {
  if (missing(frame) || missing(vsum)) stop("syntax error")
  if (class(frame) != "AEP.frame") stop("the parameter 'frame' must be of type 'AEP.frame'")
  if (is.numeric(vsum) == FALSE) stop("the parameter 'vsum' must be of type 'numeric'")  
  prob <- 0
  nb_simplex <- 0
  spx <- new("AEP.HCSimplex", array(0, c(frame@nb_var)), vsum, 1, 1)
  ListOfSimplexes <- c(spx)
  while (prob < 1 & length(ListOfSimplexes) > 0 & nb_simplex < frame@iteration) {    
    spx <- ListOfSimplexes[[1]]    
    if (length(ListOfSimplexes) > 1) ListOfSimplexes <- ListOfSimplexes[2:length(ListOfSimplexes)]
    else ListOfSimplexes <- NULL
    hc_full <- new("AEP.HCSimplex", spx@b, spx@h, spx@s, spx@iteration)
    prob_full <- VHmeasure(frame, hc_full)		
    nb_simplex <- nb_simplex + 1
    if (prob_full > 0) {
      hc <- new("AEP.HCSimplex", spx@b, frame@alpha * spx@h, spx@s, spx@iteration)
      prob <- prob + hc@s * VHmeasure(frame, hc)			
      nb_simplex <- nb_simplex + 1
      new.it <- spx@iteration + 1
      # Check how many new simplex can be generated without exceeding the limit frame@iteration      
      num_of_simplex <- frame@iteration - (nb_simplex + length(ListOfSimplexes))
      if (num_of_simplex > 0) {           
        for (j in 2:frame@N) {
          new.b <- spx@b + (frame@alpha * spx@h * frame@binary_matrix[j,])
          new.h <- (1 - (frame@alpha * frame@sum_ones[j])) * spx@h
          new.s <- spx@s * frame@m[j]				
          ListOfSimplexes <- c(ListOfSimplexes, new("AEP.HCSimplex", new.b, new.h, new.s, new.it))			
        }			  
      }			
    }    
  }      
  new("AEP.result", prob, vsum, nb_simplex)
} 

AEP.d <- function(data, alpha, iteration=1, vsum) {
  if (missing(data) || missing(vsum)) stop("syntax error")
  if (is.data.frame(data) == FALSE) stop("the parameter 'data' must be of type 'data.frame'")    
  if (is.numeric(vsum) == FALSE) stop("the parameter 'vsum' must be of type 'numeric'")        
  if (missing(alpha)) alpha <- 2 / (dim(data)[2] + 1)
  else if (is.numeric(alpha) == FALSE) stop("the parameter 'alpha' must be of type 'numeric'")      
  if (alpha < 1 / dim(data)[2] || alpha >= 1) stop("the parameter 'alpha' must be in [1/d, 1)")       
  if (is.numeric(iteration) == FALSE) stop("the parameter 'iteration' must be of type 'numeric'")    
  if (iteration < 1) stop("the parameter 'iteration' must be >= 1")	  
  frame <- as.aep.frame(data, alpha, iteration) 
  resolveAEP(frame, vsum)
}

AEP.f <- function(frame, vsum) {    
  if (missing(frame) || missing(vsum)) stop("syntax error")
  if (class(frame) != "AEP.frame") stop("the parameter 'frame' must be of type 'AEP.frame'")
  if (is.numeric(vsum) == FALSE) stop("the parameter 'vsum' must be of type 'numeric'")
  resolveAEP(frame, vsum)
}


# AEP.result -- AEP.result -- AEP.result -- AEP.result -- AEP.result -- AEP.result -- AEP.result -- AEP.result -- AEP.result -- 
# AEP.result -- AEP.result -- AEP.result -- AEP.result -- AEP.result -- AEP.result -- AEP.result -- AEP.result -- AEP.result -- 

setClass("AEP.result", representation(
  prob = "numeric",
  vsum = "numeric",
  nb_simplex = "numeric"
)
)
setMethod("initialize", "AEP.result", function(.Object, prob, vsum, nb_simplex) {
  .Object@prob <- prob
  .Object@vsum <- vsum
  .Object@nb_simplex <- nb_simplex  
  .Object
})
setMethod(f = "show", signature = "AEP.result", definition = function(object) {
  cat("P[ sum <= ", object@vsum ,"] = ", object@prob, " ( x ", object@nb_simplex , ")\n")
})

# AEP.HCSimplex -- AEP.HCSimplex -- AEP.HCSimplex -- AEP.HCSimplex -- AEP.HCSimplex -- AEP.HCSimplex -- AEP.HCSimplex -- AEP.HCSimplex -- 
# AEP.HCSimplex -- AEP.HCSimplex -- AEP.HCSimplex -- AEP.HCSimplex -- AEP.HCSimplex -- AEP.HCSimplex -- AEP.HCSimplex -- AEP.HCSimplex -- 

setClass("AEP.HCSimplex", representation(
  b = "array",
  h = "numeric",
  s = "numeric",
  iteration = "numeric"
)
)
setMethod("initialize", "AEP.HCSimplex", function(.Object, or, lim, multiplier, it) {
  .Object@b <- or
  .Object@h <- lim
  .Object@s <- multiplier
  .Object@iteration <- it
  .Object
})

# TEST -- TEST -- TEST -- TEST -- TEST -- TEST -- TEST -- TEST -- TEST -- TEST -- 
# TEST -- TEST -- TEST -- TEST -- TEST -- TEST -- TEST -- TEST -- TEST -- TEST -- 
# TEST -- TEST -- TEST -- TEST -- TEST -- TEST -- TEST -- TEST -- TEST -- TEST -- 

  par(ask=FALSE)
	par(mfrow=c(1,1))
	
  my.data <- read.table("analysis_traces.txt", header=TRUE, sep=" ")

  #my.data <- as.data.frame(my.data[,2:3])

  my.nb_obs <- dim(my.data)[1]
  my.nb_var <- dim(my.data)[2]
  
  my.nb_points <- 10
  my.nb_iterations <- 10
  my.epsilon <- 0.0001
  my.max.range = 2
  # mydata cannot contain any 0 !!!!!
  # otherwise the estimation is screwed because of the P[X_k  in (b_k, b_k + h]] in equation 1.3. The problem is the "(" with b_k not included.
  # try with 2 variables X_1 and X_2 and the following values:
  # (10, 0) and (0, 50)
  # For Q( (0,0) , 50 ) we get H(0,0) - H(0, 50) - H(50, 0) + H(50, 50) = 0 - 0.5 - 0.5 + 1 = 0 !
  
  my.data[my.data == 0] <- my.epsilon

	observed.exectime <- array(0, c(my.nb_obs))
	for (i in 1:my.nb_obs) observed.exectime[i] <- sum(my.data[i,])
  
	plot( x = 1:my.nb_obs, y = observed.exectime, 
		main = "Observed execution time",
		xlab = "#run",
		ylab = "Observed Exec Time",    
		type = "l", 
		col = "red")

	ECDF <- ecdf(observed.exectime)
  
  observed.ecdf <- array(0, c(my.nb_points))
  
  
  index <- 1  
  my.steps <- seq(min(observed.exectime), my.max.range * max(observed.exectime), length.out = my.nb_points)  
  for (i in my.steps) {
    observed.ecdf[index] <- ECDF(i)
    index <- index + 1
  }
  
	plot(x=my.steps, y=observed.ecdf,
		main = "Observed ECDF",
		xlab = "#run",
		ylab = "P[ExecTime > X]",
		ylim = c(0,1),
    type="l"
	)
  
  print("now processing V-AEP")
  
	
	estimated.ecdf <- array(0, c(my.nb_points))

  
  
	my.aep.frame <- as.aep.frame(data=my.data, iteration=my.nb_iterations, bandwidth=c(0.1,0.1,0.1,0.1,0.1))
  
	index <- 1
	for (i in my.steps) {
    result <- AEP.f(frame=my.aep.frame, vsum=i)
    print(result)
		estimated.ecdf[index] <- result@prob
		index <- index + 1
	}
	
  for (i in seq(length(estimated.ecdf), 1, -1)) {
    if (estimated.ecdf[i] > 1) estimated.ecdf[i] <- 1
    if (i < length(estimated.ecdf) & estimated.ecdf[i] > estimated.ecdf[i+1]) estimated.ecdf[i] <- estimated.ecdf[i+1]
  }
  
	lines( x = my.steps, y = estimated.ecdf, 		
		type = "l", 
		col = "red")
  










